Oscillations of dark solitons in trapped Bose-Einstein condensates 
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We consider a one-dimensional defocusing Gross-Pitaevskii equation with a parabolic potential. 
Dark solitons oscillate near the center of the potential trap and their amplitude decays due to 
radiative losses (sound emission). We develop a systematic asymptotic multi-scale expansion method 
in the limit when the potential trap is flat. The first-order approximation predicts a uniform 
frequency of oscillations for the dark soliton of arbitrary amplitude. The second-order approximation 
predicts the nonlinear growth rate of the oscillation amplitude, which results in decay of the dark 
soliton. The results are compared with the previous publications and numerical computations. 



I. INTRODUCTION 

Dark matter-wave solitons in atomic Bose-Einstein 
condensates (BECs) have many similarities with dark op- 
tical solitons in defocusing nonlinear media Q, Q • Both 
entities are fundamental nonlinear excitations of the de- 
focusing nonlinear Schrodinger (NLS) equation, describ- 
ing the evolution of electric field envelopes in the con- 
text of optics, or of the order parameter (the condensate 
mean- field wavefunction) in the context of atomic BECs. 
Dark matter-wave solitons were observed in a series of 
experiments carried out with BECs confined in external 
parabolic magnetic trapping potentials 0,0, 0|, and have 
inspired subsequent investigations of their stability and 
dynamics. In particular, if a dark (black) soliton is ini- 
tially placed exactly at the center of the magnetic trap, it 
remains standing, while if it is misplaced, it starts oscil- 
lating near the center of the trap. In this respect, there 
has been a recent interest in theoretical studies concern- 
ing the frequency and amplitude of these oscillations, and 
especially on their dependence on the velocity and am- 
plitude parameters of dark solitons. 

Preliminary studies of small-amplitude oscillations of 
dark solitons were reported in [6J, [?J by using the col- 
lective coordinate approach and assumptions of the soli- 
ton's adiabatic dynamics. The frequency of oscillations 
obtained in ||| Q was found to mismatch the correct fre- 
quency as explained in |sl Wi . The main reason for 
discrepancy of results of 0, is the use of the center-of- 
mass quantity (the Ehrenfest Theorem), as the integral 
uantity can change due to radiation from dark soliton 
Another integral quantity (the renormalized power 
invariant) was used in 01 (based on results of H3)) but 
the perturbation theory led to many correction terms in 
the main equations and relatively poor a gree ment with 
numerical simulations (see, e.g., Fig. 1 in 10]). 

Another version of the perturbation theory for dark 
solitons can be developed with the use of the renormal- 
ized momentum [l2| and can be adapted to incorporate 
radiative effects on the adiabatic dynamics of dark soli- 
tons |13| . This version of the perturbation theory relies 
on the completeness of eigenfunctions of the linearization 
problem established in 0, Ell ■ 



A recent application of 



the perturbation theory, also accounting for radiative ef- 
fects, has been reported in 16] for dark solitons in the 
presence of linear gain and two-photon absorption. The 
frequency of small-amplitude oscillations of dark solitons 
near the center of the trapping potential (and additional 
localized impurities) can be found from an integration of 
the renormalized momentum |9j in a much simpler form, 
as compared to the one presented in jlQj . 

Other integral invariants have also been used to study 
the adiabatic dynamics of dark solitons, e.g. the 
boundary-layer integral for the corresponding hydrody- 
namic equations $] and the energy (Hamiltonian) of the 
dark soliton 0] . Oscillations of dark solitons were also 
studied in the shallow soliton limit with the use of the 
Korteweg-de Vries (KdV) approximation [18(. Surpris- 
ingly enough, it was found that the frequency of oscil- 
lations of dark solitons does not depend on the soliton 
amplitude 17], or, in other words, the frequency of os- 
cillations remained the same in the limits of black and 
shallow solitons [ToL Il8| . 

Recent numerical studies of oscillations of dark soli- 
tons in parabolic and other external traps were reported 
in |2(J, |2l| . It was found that the dark solitons emit 
radiation (in the form of sound waves) due to oscillations. 
If radiation escapes the trap (as in the case of a tight dim- 
ple trap 0], or in an optical lattice potential |2l|). the 
energy (momentum) of the dark soliton changes, result- 
ing in the growth of the oscillation amplitude and the 
decay of the soliton amplitude. A phenomenological ex- 
planation of the radiative decay of the soliton energy and 
the quadratic dependence of the energy decay rate on the 
soliton acceleration was proposed in [!9ll2Q| . based on the 
earlier analysis of 0]. However, the authors of con- 
sider instability-induced dynamics of dark solitons in a 
homogeneous system, which is a different problem from 
dynamics of dark solitons in a trapped condensate. 

Further variants of the problem include oscillations of 
ring dark solitons and vortex necklaces in two dimensions 
(i.e., for "pancake" BECs) , dynamics of shallow dark 
solitons in a gas of hard core bosons (the so-called Tonks- 
Girardeau gas which, in mean-field picture, is described 
by a quintic nonlinear Schrodinger equation) para- 
metric driving of dark solitons by periodically modulated 
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Gaussian paddles |24j , and scattering of a dark soliton on 
a finite size obstacle [25|. It should be mentioned that 
the technique presented in [24| may pave the way for ob- 
serving long-lived oscillating dark matter-wave solitons 
in future BEC experiments. 

The present paper deals with the oscillations of a 
dark soliton in a BEC confined in a parabolic trap, as 
well as the inhomogeneity-induced emission of radiation. 
Our analytical approach relies on a systematic asymp- 
totic multi-scale expansion method, based on ideas of 
[l2l PHI, as well as the perturbation theory for dark soli- 
tons [l4l Il5j . We prove that the main equation of mo- 
tion for adiabatic dynamics of dark solitons of any am- 
plitude is given by the harmonic oscillator equation with 
a constant frequency. Additionally, we account for radia- 
tion escaping the dark soliton and compute the nonlinear 
growth rate of the oscillation amplitude versus the am- 
plitude of a dark soliton. Our results give a systematic 
basis for analysis of the dynamics of dark solitons in other 
systems. The analytical results are found to be very sim- 
ilar to the ones reported in recent studies 0, 0] , except 
that a regular asymptotic technique, based on the ex- 
plicit small parameter of the problem, replaces previous 
qualitative estimates based on numerical observations. 
Furthermore, we show why a formal application of the 
perturbation theory fails to incorporate the correct de- 
pendence of frequency of dark soliton oscillations. 

The paper is organized as follows. Section 2 formulates 
the problem and reports the main results. Section 3 de- 
scribes the asymptotic limit for the ground state of the 
parabolic potential. Section 4 gives a transformation of 
the Gross-Pitaevskii (GP) equation to the regularly per- 
turbed nonlinear Schrodinger (NLS) equation. Section 5 
contains the analysis of the perturbed NLS equation up 
to the first-order and second-order corrections. Section 
6 reviews the formal application of perturbation theory 
to the GP equation and the previous results. Section 7 
concludes the paper. 

II. MODEL AND MAIN RESULTS 

At low temperatures, the dynamics of a repul- 
sive quasi-one-dimensional BEC, oriented along the x- 
axis, can be descibed by the following effective one- 
dimensional (ID) GP equation (see, e.g., j2(|) 

h 2 

ifK/j t = -—tp xx + V(x)^ + g\i;\ 2 iP, (1) 
Zm 

where subscripts denote partial derivatives, ij){x,t) is the 
mean-field BEC wavefunction, m is the atomic mass, and 
the nonlinearity coefficient g (accounting for the inter- 
atomic interactions) has an effective ID form, namely 
g = 2hauj±, where a is the s-wave scattering length and 
lu±_ is the transverse-confinement frequency. Addition- 
ally, the external potential V(x) is assumed to be the 
usual harmonic trap, i.e., V{x) — rau) 2 x 2 /2, where u> x is 
the confining frequency in the axial direction. 



To reduce the original GP equation to a dimension- 
less form, x is scaled in units of the fluid healing length 
£ = h/^/nogm (which also characterizes the width of the 
dark soliton), t in units of £/c (where c = W nog/m is the 
Bogoliubov speed of sound), the atomic density n = \ip\ 2 
is rescaled by the peak density no, and energy is mea- 
sured in units of the chemical potential of the system 
fj, = griQ. This way, the following normalized GP equa- 
tion is readily obtained, 

iik — -~u xx + e 2 x 2 u + \u\ 2 u, (2) 

where the parameter e = (2y/2ano)~ 1 (uj x /uj±) deter- 
mines the magnetic trap strength and u(x,t) 6 C. Let 
us assume realistic experimental parameters for a quasi- 
1D repulsive condensate containing N ~ 10 3 -10 4 atoms 
and with peak atomic density no ~ 10 s m _1 . Then, as 
the scattering length a is of order of a nanometer (e.g., 
a = 5.8 nm or a = 2.7 nm for a 87 Rb or 23 Na conden- 
sate), and the ratio of the confining frequencies (for such 
a quasi-lD setting) is ui x /oj± ~ 1/200, it turns out that 
the magnetic trap strength e is typically O(10~ 2 ). Thus, 
e is a natural small parameter of the problem. 

When e = 0, the defocusing GP equation J5J has the 
exact solution for the dark soliton: 

u ds (xA) = [fctanh(A;(x - vt - s )) + iv] e ~ lt+ie ° , (3) 

where k = \/l — v 2 < 1 is the amplitude of the dark soli- 
ton (with respect to the continuous- wave background), 
\v\ < 1 is the velocity parameter, and (sq,9q) G M 2 are 
arbitrary parameters of the position and phase. Since 

|wds| 2 = 1 — fc 2 scch 2 (/c(a; — vt — s)), 

it is clear that the continuous-wave background for the 
dark soliton (or the dimensionless chemical potential fio) 
is normalized by one, such that lim|a.|_ >cx) |itds| 2 = 1- 
When k — > 1 and |k| — * 0, the dark soliton approaches the 
limit of a standing topological soliton (called the black 
soliton). When k — > and |u| — > 1, the dark soliton ap- 
proaches the limit of a small-amplitude shallow soliton 
(which satisfies the KdV approximation 18]). 

It should be noticed that, in physical terms, the 
choice no — 1 actually sets the number of atoms N 
of the condensate. In the framework of the Thomas- 
Fermi approximation |27|. it can be found that N — 

(4\/2/3)(^no/0)/io /2 , and, thus, for n « 10 s m -1 , 
£ ~ 0.1-1 microns (which are realistic value of the heal- 
ing length for quasi-lD Rb BECs) and Q - 10~ 2 , the 
choice /iQ — 1 leads to a number of atoms of the order of 
N ~ 10 3 -10 4 . 

When e / 0, but the nonlinearity is crossed out by the 
linearization, the parabolic potential of the GP equation 
has the ground state solution: 

/ ±\ ( ^x 2 + iet\ 
u gs (x, t) = u exp -= — , (4) 
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where uq £ C is an arbitrary parameter of the ground 
state amplitude. When the nonlinear term of the defo- 
cusing GP equation (J2J) is taken into account, the linear 
mode (@J) generates a family of ground state solutions by 
means of a standard local bifurcation 1281. 



(5) 



where /i e e R is the normalized chemical potential and 
6>o £ K is an arbitrary phase. Due to the scaling invari- 
ance of the GP equation (J2J , the amplitude of the ground 
state U e (x) £ R can be uniquely normalized by one, such 
that \U C (0)\ 2 = 1. 

The first excited state of the parabolic potential bifur- 
cates from the linear solution: 



u les (x, t) = xexp 



ex 



3iet 



V2 



(6) 



by means of the same local bifurcation 28]. The first 
excited state corresponds to a static bound state between 
the dark soliton JSJ placed at the center x = of the 
nonlinear ground state (JSJ. The solution for the static 
bound state exists for any e ^ but tells nothing about 
dynamics of the dark soliton placed near the center of 
the ground state. 

Dynamics of dark soliton is considered in this paper. 
We show that the dark soliton © undertakes adiabatic 
dynamics in the limit e — > 0, such that the parameter 
(sq + vt) = s(T)/e of position of the dark soliton © 
becomes a function of slow time T — et, while the ve- 
locity parameter v is defined as v(T) — s. The adia- 
batic dynamics leads to generation of radiative waves, 
which escape the dark soliton but become trapped by 
the parabolic potential. In the decomposition of the solu- 
tion u(x,t) into two (inner and outer) asymptotic scales, 
the leading-order radiative effects are taken into account 
when parameter 0q = S{T) of complex phase of the dark 
soliton © depends also on T = et and the first-order cor- 
rections to the dark soliton © grow linearly in x. When 
reflections from the trapping potential are neglected, the 
extended dynamical equation for the position s(T) of the 
dark soliton © takes the form: 



cs 



2^/(1 - s 2 )Vl - s 2 - s 2 



0(e 2 



(J) 



in the domain (s,s) £ T> , where T> is the unit disk: 

V = {{8,a) £ M 2 : s 2 +s 2 <l}. (8) 

The left-hand-side of the dynamical equation Q repre- 
sents the leading-order adiabatic dynamics of the dark 
soliton oscillating on the ground state of the trapping 
potential. The right-hand-side represents the leading- 
order radiative effects (sound emission), when reflections 
of radiation from the parabolic potential are neglected. 
Within the truncation error of 0(e 2 ), the ground state 
iJSJ can be approximated by the Thomas-Fermi (TF) ap- 
proximation U e {x) = yl — e 2 x 2 on the scale \x\ < e^ 1 . 



The only equilibrium point of the dynamical equa- 
tion (2J is (0,0) and it corresponds to the static bound 
state bifurcating from the first excited state ©■ The 
leading-order part of the dynamical equation (JJJ de- 
scribes a harmonic oscillator with the obvious solution: 
s(T) = sq cos(T+8n). This result is well-known from ear- 
lier papers y| M> HH 1 where it was derived in the limit of 
black soliton, when k — > 1 and \v\ — > 0. In the derivation 
presented herein, we have not used the assumption on the 
initial position s(0) and speed s(0) of the dark soliton, 
and therefore, the approximation of the harmonic oscilla- 
tor Q) remains valid for larger values of (s, s) inside the 
unit disk ©. 

We note that the velocity-dependent correction to the 
frequency of the harmonic oscillator Q was obtained in 
an earlier study (see Eq. (36) in [10|), but it is not con- 
firmed within our analysis. On the other hand, our re- 
sults confirm the results of the shallow soliton approxi- 
mation |l8j which establishes the same frequency of os- 
cillations as in the black soliton limit. In addition, the 
uniform frequency of oscillations for dark solitons of all 
amplitudes and velocities was recently reported by means 
of the energy (Hamiltonian) computations |17| . 

Let E be the energy of the harmonic oscillator: 



E=-{s 2 + s 2 ) 



(9) 



The energy increases in time due to the first-order cor- 
rection terms of the main equation (|7|): 



E 



vT 



= + O(e 2 )>0, (10) 



where (s, s) £ T>q. Due to the energy pumping l |10|l . the 
amplitude of the harmonic oscillator increases in time. 
When the oscillation amplitude grows, the dark soliton 
© shifts from the black soliton limit k — ► 1, \v\ — > to 
the shallow soliton limit k — ► 0, \v\ — ► 1. By the Poincare- 
Bendixon Theorem, the limit cycle does not exist in the 
unit disk JSJ and all orbits approach the boundary of the 
disk, where the main equation |Q becomes invalid. The 
growth rate of the oscillation amplitude is nonlinear in 
general and depends on (s, s). 

In the limit s 2 + s 2 — » 0, the main equations Q and 
(|10|l can be simplified. First, the energy of the dark soli- 
ton oscillations accelerates by the squared law E = es 2 /2, 
which is postulated in [Tiil ] and confirmed in numerical 
computations where the dark solitons oscillated in a tight 
dimple trap (see Fig. 2 in 0]). Second, the nonlinear 
equation J7J) is linearized as follows: 

s + s-|s = 0(e 2 ,s 3 ), 

which shows that the center point (0, 0) becomes an un- 
stable spiral point on the plane (s, s) as < e -C 1, with 
the leading-order solution: s(T) — soe tT ^ 4 cos(T + S a ). 
Therefore, due to radiative losses, the amplitude of oscil- 
lations of dark solitons increases while its own amplitude 
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decreases. This main result of the asymptotic analysis 
was also confirmed by numerical simulations where the 
dark solitons oscillated between two Gaussian humps (see 
Fig. 9 in Hj). 

Figure ^ shows a spatio-temporal contour plot of the 
reduced density (the ground state density subtracted 
from the actual density) for a one-dimensional BEC 
confined in a parabolic trap with normalized strength 
e = 0.05 (the TF radius is equal to 20). The figures 
were obtained by numerical integration of the GP equa- 
tion J21 . The white areas correspond to a dark soliton, 
initially placed at the trap center (i.e., s(0) = 0) with an 
initial velocity v(0) = 0.1 (bottom panel) and v(Q) = 0.5 
(top panel). The plot is compared to the two versions 
of the main equation 10. The solid lines correspond to 
harmonic oscillations (within the adiabatic soliton dy- 
namics), while the dashed ones correspond to the grow- 
ing oscillations (within the inhomogeneity-induced sound 
emission) . It is seen from the figures that the dashed lines 
approximate better the actual dark soliton motion within 
the first period of oscillations for < t < < T osc , where 
T osc — 2n/e. For instance, a better agreement is achieved 
at the first turning points where, due to stronger sound 
emission, a slight increase of the amplitude of soliton os- 
cillations is readily observed. 

Nevertheless, for longer times, the anti-damped ap- 
proximation of the main equation (0 becomes irrelevant 
due to the fact that the radiation cannot escape the trap 
and is reflected back to the dark soliton. As a result, 
the soliton continuously interacts with the reflected radi- 
ation so that, on average, the soliton reabsorbs the radi- 
ation it emits on the contrast to numerical computations 
in 0, [2(j. We note that the sound emission is much 
weaker for the deeper soliton (with v = 0.1) and conse- 
quently the analytical result pertaining to the adiabatic 
approximation of the soliton motion is much closer to the 
result obtained by the numerical simulation. 

The developed method allows us to incorporate mul- 
tiple sound reflections and their recombination to the 
asymptotic equations for the dark soliton. If this is done, 
the main equation is extended by an additional term 
due to multiple reflections beyond the initial time inter- 
val < t < £*. This complication is however beyond the 
scopes of the present paper. 

III. GROUND STATE OF THE GP EQUATION 

We start with analysis of the ground state solution of 
the GP equation in the form 0, where (U £ {x),n e ) 
is a real- valued pair of eigenfunctions and eigenvalues of 
the nonlinear boundary- value problem: 

-U" - eVt/ - U 3 + fiU = 0, x £ R+, (11) 

subject to the normalized boundary conditions: 

U e (0) = 1, U' e (Q) = 0, lim U e (x) = 0. (12) 
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FIG. 1: Spatio-temporal evolution of the reduced condensate 
density (the ground state density minus the actual density) for 
the GP equation with e = 0.05, s(0) = 0, and v(0) = 0.1 
(top two panels), v(0) = 0.5 (bottom two panels). The solid 
lines correspond to the solutions of the harmonic (adiabatic) 
approximation and the dashed lines correspond to the anti- 
damped approximation (pertaining to sound emission). The 
dotted line on the plot of \u{x, 100) | 2 shows the parabolic 
trapping potential. 



We are interested in existence of the symmetric ground 
state U e (x) > 0, U e (— x) = U e (x) on x e R for small 
values of e 2 . We recall two facts for the boundary- value 
problem (flip with a given value of e > (see |13|): (i) the 
local bifurcation from the linear ground state 10 occurs 
for fi > /io( e )j where /io = c/V2, and (ii) the nonlinear 
ground state exists as a one-parameter smooth family, 
parameterized by [i or, equivalently, by t/(0). Moreover, 
(7(0) is an increasing function of fi > /J-o(e) for a fixed 
value of e. Therefore, for a given value of e > 0, there 
exists a unique value of [i (called fi e ), which corresponds 
to the normalization U e (0) = 1 and the solution pair 
(U t (x), /i e ) is a smooth function of e. The solution of 
the boundary- value problem of ljll|l - H12|l can be approx- 
imated numerically by means of a contraction mapping 
method. Figure 2 shows the dependence of /j, e versus e, 



5 





FIG. 2: Ground state solution of the boundary- value problem 
l|lllL Top panel: The dependence of /j, e versus e for fixed 
U e (0) = 1. Middle panels: Profiles of the function U(x) for 
two different values of the trap strength: e — 0.025 (left) and 
e — 0.25 (right). Bottom panel: The dependence of £7(0) 
versus ft for fixed e = 0.05. 



the profile U e (x) for different values of e, and the depen- 
dence of U(0) versus /x, for fixed e = 0.05. 

The solution of the ODE Hll|) with initial values 
U(0) = 1 and t/'(0) = can be constructed in the power 
scries form: 



U(x) 



where coefficients {a/c}£° 
the first two terms are: 



E 

k=l 



ClkX 



2 A- 



(13) 



can be found recursively, e.g. 



ai = 1 - fi, 

a 2 = i (e 2 + (1 - a*)(3 - M )) . 

The parameter /i = /i e is defined from the decay condi- 
tion at infinity, where U = U e (x) and lim^^oo U e (x) = 0. 
However, the decay condition is computationally ineffi- 
cient for approximations of the dependence of /i e versus e. 
A different (WKB) method was used for approximations 
of the solution U e (x) and the dependence [i t in nam. 



The advantage of the WKB method over the power se- 
ries <|13|) is that the dependence /i e is determined from the 
condition at x = rather than from the decay condition 
as x — > oo. 

In the WKB method, the solution of the ODE l(TT|) 

is represented in the form U{x) 
satisfies the equivalent problem: 



Q(x) = fi 



2 2 i 

e x + 



\/Q{x), where Q(x 
)"(x) Q'\x) 



4Q(x) 



(14) 



In the limit of small e, the solution of the ODE 114|) can 
be thought in the form of a WKB asymptotic series: 



Q = M - X" 



k(X), 



(15) 



where X = ex, \X\ < ^/jl, and the set {Qk{X)}^° =1 can 
be found recursively, e.g. the first two terms are: 

/' 



Qi 



2(M-^ 2 ) 5 



and 



}'{{X) , XQ' l {X)+Q 1 {X) X 2 Q,(X) 



4( M -X2) 



2(p~X 2 ) 2 



in-x*r 



The symmetry in X implies that the condition Q'(0) = 
is satisfied. The normalization condition Q(0) = 1 
defines the power approximation of the dependence /i = 
fx e versus e: 



e 



such that 



e 2 e 4 
1 + T + 2 



0(e 6 ) = l, 



0(e 6 ). 



(16) 



(17) 



Using the leading-order approximation for we find 
an asymptotic approximation for U e (x) from the WKB 
asymptotic series (|15fl : 



U e (x) = v 7 ! -e 2 x 2 + € 2 U(ex,e), 



(18) 



where U(ex,e) is the remainder term. The leading-order 
approximation in IjlSB is referred to as the Thomas- 
Fermi approximation |27l |. The WKB asymptotic series 
(|15f) gives a bounded approximation of the solution pair 
(U t {x), /z £ ) for \ex\ < ^fJT t but it becomes invalid at and 
beyond the turning points at \ex\ >yfjl^ |29| . Therefore, 
the Thomas-Fermi approximation l|18|) is valid only for 
lea; I < 1 in the limit of small e. 



IV. TRANSFORMATION OF THE GP 
EQUATION 

We proceed with analysis of dynamics of the dark soli- 
ton on the ground state of the parabolic potential. This 
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problem can be studied after the equivalent transforma- 
tion of the GP equation J2J) : 



u(x,t) = U t (x)w(x,t)e~ Hl ' t , 



(19) 



where (U e (x), jjL c ) is the solution pair of the boundary- 
value problem (|ll[l - i|12[) and w(x,t) is a new variable. 
The function w(x, t) satisfies the PDE problem: 

iwt + \w xx + U 2 (x)(l - \w\ 2 )w = ~ J ffj^Wx- (20) 

We will use the fact that U e (x) = U e (X), X — ex in 
the asymptotic region \ex\ < 1 as e — > 0. Assuming 
that the dynamics of the dark soliton occurs inside the 
asymptotic region \ex\ < 1, we introduce the traveling 
wave coordinate for the dark soliton: 



r\ = x 



s(T) 



T = et, 



(21) 



where the dependence s(T) is to be determined, and 
expand the function U e (X) in the Taylor series near 
X = s(T): 

U e (X) = U e (s(T)+er}) 

= U e (s) + e V U' £ (s) + \eWU':{ S ) + 0((er,) 3 ). (22) 

As a result, the PDE problem (|20() takes the form of the 
perturbed equation: 



lW t - IVWr, + ^W nn + U 2 (s)(l - \w\ 2 )lV 



1 

2 



w r) + 2U e (s)U' t (s)i 1 (\ ~ \w\ 2 )w 



2 {U':{ S )U t { S )-U' 2 {s)) 



U?(s) 



- 2 rf (UMU'Jis) + U' 2 (s)) (1 - \w\ 2 )w + 0(e 3 ), 

where v(T) — s is the speed of the dark soliton. The 
perturbed equation is simplified with the scaling trans- 
formation: 



z = vU e (s(T)), v(T) = v(T)U e (s(T)), 



(23) 



such that the wave function w{z, t) satisfies the perturbed 
defocusing NLS equation: 



iw t + U 2 {s) 



-ivw z + l;W zz + (1 - \w\ 2 )w 



+R(w,w) = 0, (24) 



where R = eRi + e 2 i?2 + 0(e 3 ) are perturbation terms, 
with the first two of them being: 

Ri = U' t {s){ivzw z +w z + 2z{l-\w\ 2 )w) 
(U?(a)U e (a)-U?( 3 )) 

Ri = TTTT^-, ZW Z 



U 2 {s) 

{U t { S )U':(s) + U' 2 {s)) 
U 2 {s) 



z 2 {\ - \w\ 2 )w. 



Since U e (s) > for the ground state of the parabolic 
potential, the perturbation terms of the equation l|24JI are 
regular for any s£K. However, since the representation 
U e (x) = U e (X) is valid only for \X\ < 1 in the limit of 
small e, we consider the perturbed NLS equation l|24() in 
the region where \s\ < 1. The leading-order part of the 
defocusing NLS equation (|24|l (when R(w, w) — 0) has 
the exact solution for the dark soliton: 



w{z, t) = wq(z) — k tanh(/«z) + iv, 



(25) 



where k = \f\ — v 2 and v(T) — vq is a constant speed. 
The time evolution of the dark soliton (|25l) in the case 
R(w, w) y^z is studied with a reg ular perturbation theory 
for dark solitons LL3, [13, llj, LLal . 



V. ASYMPTOTIC MULTI-SCALE EXPANSION 
METHOD 

The solution to the perturbed NLS equation (PU) can 
be obtained in the form of an asymptotic multi-scale ex- 
pansion series for the perturbed dark soliton i12ol) 13]: 



w 



(z, t) = [to + ewi + e 2 w 2 + 0(e 3 )] e 



i0 



(26) 



where T = (et, e 2 t, ...), the function wq = wq(z; T) is the 
dark soliton 125|) . while the functions w\ — wi(z,t;T) 
and ui2 — W2(z,t;T) solve the inhomogeneous linear 
problems: 

i<9t(73Wi + U 2 (s)TLwi = 9wq — idrcr^o 

-RiK,wo) (27) 

and 

idt<T3W2 + U 2 (s)7iw2 = 0Wi — idxc^Wi 
-N 2 (w 1 ,w 1 ) - DRi(w ,u)o)wi - R 2 (w ,^o)- (28) 

Parameters s(T) and 9(T) are to be determined from 
solutions of the inhomogeneous problems, while v(T) = 
s/U e (s). In the problems and we have intro- 

duced the notations: and Rfc for vectors (wk,u>k) T , 
k = 0, 1, 2 and (Rk,Rk) T , k = 1,2, DUi for the Jacobian 
of Ri , N2 for the quadratic terms from the left-hand-side 
of the unperturbed NLS equation. The self-adjoint lin- 
earization operator is 



H = —ivuzdz 



l - 



2\w \ 



2\w \ 



where <tq = diag(l, 1) and 03 = diag(l, — 1). Analysis 
of the first-order problem 127|) predicts a leading-order 
equation for s(T), which characterizes oscillations of dark 
solitons near the center of the parabolic potential. The 
first-order problem also describes the leading-order radia- 
tion from the dark soliton to infinity, related to the equa- 
tion for 9(T). Analysis of the second-order problem l|28|) 
predicts a first-order correction to the equation for s(T), 
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which is induced by the leading-order radiation. Due to 
radiation, the oscillation amplitude grows in time so that 
the amplitude of the dark soliton decreases. Derivation 
of all these results is divided into several technical sub- 
sections. 

A. Linearization operator 



The eigenvectors JHTJ and © are not in L 2 (R, C 2 ). This 
fact implies (see 13] ) that adiabatic dynamics of the dark 
soliton (the decaying component) induces radiative waves 
of the continuous spectrum (the bounded non-decaying 
component) already at the first order of the perturbation 
theory. The coupling between the decaying and non- 
decaying components is computed from the dynamical 
equations on parameters s(T) and 9(T). 



The self-adjoint operator TL is defined on complete 
Hilbert space H l {R, C 2 ) C L 2 (K,C 2 ). It has a non- 
empty kernel: 



Hw' (z) = 0, 



(29) 



which is related to translational invariance of the unper- 
turbed NLS equation in z. Furthermore, the operator TL 
has two branches of the continuous spectrum: 



S (W) = (-oo,-2] U (-oo.O], 



(30) 



such that the second branch intersects the kernel. The 
only bounded non-decaying eigenvector for the zero 
eigenvalue, which belongs to the continuous spectrum 
(|3^|l , is related to the gauge invariance of the NLS equa- 
tion: 



TL (ia 3 w ) = 0. 



(31) 



The homogeneous equation TLw = supports the decay- 
ing solution (|29[) . the bounded solution l|31[l . a linearly 
growing solution, and an exponentially growing solution. 
The linearly growing solution can be found explicitly ^3 : 



TL (icr 3 zw - <9„w + T^A w o^ = 0. (32) 

Here and henceforth, it is convenient to consider wo(z) as 
a function of two independent parameters n and v. The 
relation k = yl — v 2 is used after evaluating the partial 
derivatives of wq(z) in k and v. 

The linearization operator a 3 TL has the kernel (|29|l and 
the generalized kernel: 



c 3 W ^„w <9 K w J = iw„(z) 



(33) 



The spectrum of a 3 TL includes two branches of the con- 
tinuous spectrum: 



a css {a 3 H) =1U 



(34) 



The spectrum of a%H consisting of the kernel (|29|l . 
the generalized kernel (|33|l and the two branches of 
the continuous spectrum (|34J) is complete in i? 1 (IR, C 2 ) 
[T4I llfij. Nevertheless, since the kernel of TL has also a 
bounded non-decaying eigenvector l|31|l . we construct a 
bounded non-decaying eigenvector from the nonhomoge- 
neous problem: 



03H [ ^r^ w o ) = i (icr 3 w ) . 



(35) 



B. The leading-order frequency of oscillations 

The linear inhomogeneous problem l|27|l leads to a sec- 
ular growth of wi(z,t;T) in t unless the right-hand-side 
is orthogonal to the kernel of TL. The orthogonality con- 
dition defines a nonlinear equation on parameters of the 
dark soliton 125fl : 



1 



sech 2 (K,z)lm(dTWo)d(nz) 



sech 2 («:z)Re(i?i(wo, wq))cI(kz). (36) 



By computing the integrals directly, we obtain the system 
of dynamical equations: 



v={l- v 2 )U' t {s), 



(37) 



where the last equation is due to the relation l|23|l be- 
tween v(T) and v(T) — s. Closing the system of dynam- 
ical equations, we find the governing equation for the 
position of the dark soliton: 



s + V'(s) = 0, V (s) = -(1-U?(s)) 



(38) 



The governing equation (|38J) is equivalent to the Hamil- 
tonian system of a particle moving in a potential field, 
where V(s) stands for the effective potential energy of 
the particle. 

Using the WKB asymptotic series (|15|) for U 2 (x) and 
the power approximation (|17Jl for /x £ , we find the power 
approximation for the potential function V(s): 



e 2 s 2 (2- 



y(s) = y + 4(i - s 2 ) 2 



0(e 4 



(39) 



where |s| < 1. Since the system (|38|l is valid at the 
leading order of the asymptotic series, the dark soliton 
oscillates as a harmonic oscillator in the limit e — * 0: 



s = 0, 



(40) 



which is equivalent to the left-hand-side of the main 
equation J7J. The equation l|4UI) for harmonic oscilla- 
tor is valid in the unit disk (JSJ. Since the solutions 
s(T) — s cos(T + (5 ) represent circles of radius s on 
the phase plane (s, s), the trajectories remain inside the 
disk JSJ) whenever sq < 1. 
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We note that the 0(e 2 ) corrections of the power ap- 
proximation H39|) are beyond the Tomas-Fermi approxi- 
mation p8|l. These terms can be dropped in the main 
equation (|38[) . since the 0(e) corrections to the main 
equation l|38[l occur from the solution of the second-order 
inhomogeneous equation (|28H . 



C. The first-order radiation corrections 

After the constraint 1(36(1 is added, the linear inho- 
mogeneous equation 1(2*71) can be solved for w\(z,t;T), 
such that w\(z,t;T) is bounded in t. We adopt a stan- 
dard assumption of the asymptotic multi-scale expansion 
method that the solution w\(z,t;T) approaches the sta- 
tionary solution wi s (z;T) as t — > oo due to dispersive 
decay estimates. The stationary solution wi s (z;T) can 
be represented in the form: 



q(T) 

U?(s) 

Zvq(T) - 9(T) 



2k[/ £ 2 ( 



(izwo - d^wo) 

d K w Q +wis(z;T). 



(41) 



The first two terms in (|41|1 represent the linearly grow- 
ing ((32(1 and bounded 1(3511 eigenvectors. The last term 
wi s (z]T) is a ^-independent solution of the inhomoge- 
neous equation (|27|) which corresponds to the last two 
terms in the right-hand-side of 1(27(1 under the constraint 
(|5T|l . We do not include the decaying ("2"""|) and bounded 
(131(1 solutions of the kernel of 7i in the representation 
(|41|l as they renormalize parameters 9(T) and s(T). 

The dependence of q(T) and 9(T) is defined from the 
radiation problem, associated to the behavior of the sta- 
tionary solution w\ s (z;T) as \z\ — > oo. Let wf s (z;T) 
represent a linearly growing solution w\ s (z;T) as z — > 
±oo. Neglecting exponentially small terms in the limits 
z — ► ±oo, we obtain from l(27|) the linear inhomogeneous 
problem for wf s (z; T): 



,.±" 



ivwfg — wf s — (k ± ivYwi 



9(iv ± k) H — (k ± iv). 



(42) 



The most general linearly growing solution of the inho- 
mogeneous equation "4~2l has the form: 



1 



[(at ±b x )iz(±K + iv) + (<z 2 ± h)] , (43) 



where a\ 2 and b\ 2 satisfy two relations: 



2nb 2 = 



vb\ — 2na,2 = 



(44) 



The first two terms in the solution 1(41(1 have odd real 
parts and even imaginary parts in z, while the component 
w\ s (z,T) has the opposite symmetry in z. Matching the 



linear growing terms in (|41() and (|43|l under the relation 
(|4~""}) . we obtain that 



vq-9 
2k 



(45) 



We note that the constant terms in 141fl and l|43|l are 
equivalent to the second equation (1451) under the renor- 
malization of 9 = 9^ + e9f + 0(e ), where 9f = ±q/n. 
One more equation is needed for finding of values of a 2 
and 61. This equation can be derived from the balance 
equation for the renormalized power of the perturbed 
NLS equation J2H): 

(d t - vU 2 (s)d z ) n(w,w) + U 2 (s)8 z j(w,w) = l(w,w), 
where 

n = \w\ 2 - 1, 
j = — (ww z - w z w) , 



2 1 



I = i (wR(w,w) — wR(w,w)) . 



(46) 



Computing explicitly the integral quantities and the 
jump conditions: 



N = 
L = 



n(w, w)dz = —2k + 0(e), 

o 

l(w, w)dz = 2vKU' t (s)e + 0(e 2 ) 



and 



U 2 (s) [n(w, w)]t = 4Ka 2 e + 0(e 2 ), 
U 2 (s)[j(w,w)}t=2b 1 e + 0(e 2 ), 

we find another relation on a 2 and b\: 

Kbi - 2K 2 va 2 = vK 2 U' t (s) - vv. (47) 

Using the leading-order equation l|37|) for v, we obtain 
that 



b\ = 2kv<12 = U'(s). 

K 



(48) 



Once the parameters a±, a 2 , b\, and 6 2 in the asymptotic 
representation (|43f) are uniquely found, the stationary 
solution w s (z;T) can be defined outside the dark soliton 
up to the order of 0(e 2 ) terms: 



lim w s (z,T) = (l + eW ± (X,T))e i&± ^ T) 

z— >±oo 

= [(±K + tv) + ewt(z;T) + 0(e 2 )] e *(«o ± +^i ± +0(* a )) ) 
where X = ex, ez = U e (s(T))(X - s(T)), 



w± = ^±o 2 ) +0(£); 



U 2 (s) 

e ± (T) + e ef(T) + ez 



(gi±ji) 

U 2 (s) 



0(e 2 ), 
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and the explicit forms for and &f are not written. 
The radiation fields W ± {X, T) and Q ± (X, T) are defined 
outside the dark soliton for X ^ s(T), respectively, sub- 
ject to the boundary conditions: 



W ± 

d@ 
dX 



X=a(T) 



X=s(T) 



n(b 2 ± a 2 ) 

~ufW~ 

(oi±6i) 



+ 0(e), 



Ue(s) 



+ 0(e). 



(49) 



These conditions match the inner asymptotic expansion 
for z — 0(1) and the outer asymptotic expansion for 
X = 0(1) in the stationary solution w s (z; T) (see 



D. Radiation problem for small-amplitude waves 

We now consider the small-amplitude waves within the 
original equation 120(1 for w(x,t). Using the polar form, 
w(x, t) = R(x, t) exp (i$(x, i)), we obtain the system: 



Rt + R x ^x + ^R^xx 



uM 

U e {x) 



R&X 



R.i 
2R 



®t + ^i-^-U?(x)(l~R 2 ) 



U' e (x) Rx 

U e (x) R 



The small-amplitude long-wave solutions of the above 
system can be constructed in the asymptotic form: 

R = l + eW{X,T) + 0(e 2 ), 
<t> = 9(X,T)+0(e), 

where X = ex, T = et and we use the fact that 
U e (x) = U e {X) in the asymptotic region \ex\ < 1 as 
e -> 0. The leading-order terms W(X,T) and 0(X,T) 
solve the coupled problem: 



W T + (U e (X)V) x + 2U' e (X)V 



0. 



(U e (X)V) T + (U 2 (X)W) x = 



where 



V 



e 



X 



2U e (X) 



(50) 
(51) 



(52) 



Equivalently, the coupled problem (15011 15 1|) with the cor- 
respondence (|52|l reduces to the wave equation with a 
space-dependent speed: 



9tt- {U 2 {X)Q x ) x =0 



where 



W 



On 



2U*(X) 



(53) 



(54) 



The system (|50f) - (|51|l and the scalar equation (|53[) are 
to be solved separately for X > s(T) and X < s(T) 
subject to the boundary conditions (|49|l . In addition, 



two radiation boundary conditions must be added to the 
system (f3U|) — (I^Tl) for a unique solution. Since the small- 
amplitude waves move faster than the dark soliton (in- 
deed, s 2 < U 2 (s) = 1 — s 2 in the domain ©), the radia- 
tive waves include the right-travelling wave for X > s(T) 
and the left-travelling wave for X < s(T). Figure 
shows the upper-half {X, T)-plane, which is divided by 
the curve X = s(T) into two domains D t and D\. 

The wave equation l|53(l has two characteristics, which 
are defined by the principal part of the PDE system 
dSOJ— (EH) - Let X = £±(T) be the equations for the two 
characteristics curves, starting from a particular point 
(s(tq), Tn), where To > 0. According to the standard text 
on PDEs [3(|, we find that the functions £±(T;tq) solve 
the initial- value problem for T > tq: 

^ = ±U e (0, £±(t ;t )= S (to). (55) 

The components R± = W ± V = R±{T;t ), defined 
along the characteristics X = £±(T; r ), solve the system 
of evolution equations: 

^ = ~U' e (t + (Tm))(5R + -IL.), (56) 



dR 



ilT = --U' e (^(T;r ))(R + -5R-). (57) 

Integrating the initial- value problem (|55|l in the Thomas- 
Fermi approximation U e (X) = yl — X 2 , we find the ex- 
plicit solution: 



U(T;t ) =sin(e ±(T-7D)), 



(58) 



where £o = arcsin(s(ro)). The families of two character- 
istics intersect transversely the dark soliton curve X = 
s(T) at any point (s(t ),t ) (see Figure |3J), where the 
components R± are generated by means of the boundary 
conditions (|49|l and the radiation boundary conditions, 
which need to be added to the system (|S^|l - lf5T|) . 

Let us consider the domain D r to the right of the curve 
X = s(T). By geometry of the initial- value problem (see 
Figure |3J) or by the symmetry of the coupled problem 
(|5uTl l|57|). the characteristics ^_(T;ri) and the compo- 
nent R-(T;ti) in D r are the same as the characteristics 
£_l_(T; to) and the component i?+(T; to) in D r . When in- 
tegrating the evolution equation for R+{T; To), one needs 
to substitute the value for from its value on the 
transversally intersecting characteristics in D r . Figure 
3 shows two intersections of £+(T;to) for To < T < t\ 
with the characteristics to the point (s(r), t). Before the 
reflection, £ + (T;to) intersects with £_(T;t), such that 
i?_ = R-(T;t). After the reflection, £ + (T;to) intersects 
with £+(T; t), such that R_ = R+{T; r). At T = n, the 
matching formula gives a required radiation boundary 
condition: 



R-(tv,ti) = R+(ti,t ) 



(59) 



The initial data for i? + (To,To) are uniquely defined from 
Pty and 
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Similarly, in the domain D\, we obtain the radiation 
boundary condition: 



R+(n;n) = R-(n,T ). 



(60) 



If no incoming waves are imposed initially for T < 0, 
then we obtain that R- = in D r and R + = in D\ at 
least for < T < r* , where r* is the first intersection of 
£+(T; 0) with s(T), such that £(r*; 0) = s(r*). Therefore, 
the components V^X, T) and W A± (X, T) of the radiative 
waves to the right and left of the dark soliton X = s(T) 
are related at the moving boundary X — s(T) by the 
radiation boundary conditions: 



X=s(T) 



(61) 



X=s(T) 



The latter conditions result by virtue of (|49|l and (|52[l in 

the constraints: 



di = 2na,2, b\ = 2k& 2 . 



(62) 



Using explicit formulas 1451) and 148|) . we find unique 
expressions for q(T) and #(T): 



C/ e '( S ) 



0. 



(63) 



Relations (f Cj 1 1> — (f f>3|) are valid within the first period of 
oscillations < T < < 2tt under the condition that 
no incoming waves are generated initially. If absorbing 
boundary conditions for radiative waves are specified on 
the boundary of the TF radius (X i=s ±1), the relations 
(|6H l|63l) are extended for later times T > . Otherwise, 
additional terms occur in the relations (|61|l - l|63|) due to 
multiple reflections, and these terms are defined by the 
solution of the coupled evolution problem (|56*|) - l(5"?| along 
characteristics (|58fl with the radiation boundary condi- 
tions - . 



E. The first-order corrections to the main equation 

Radiative corrections to the main equation for s(T) are 
derived from the second-order inhomogeneous equation 
(|28|l by means of the same orthogonality condition (13611 . 
Equivalently, radiative corrections can be computed from 
the balance equation for the renormalized momentum 
M: 



(d t - vU 2 e { S )d z )p{ Wl w) + Uf(s)d z r(w,w) 
= d z m(w, w) + k(w, w), 



(64) 



where 



2 ( ww * 



ww z 



w ■ 



1 



[ww z 



2w z w z + ww zz ) 1 



2w[ 




FIG. 3: The domain of the PDE system JSDJl-(inj and the 

family of two characteristics, superposed with the dark soliton 
curve X = s(T). 



and 



m = — — (wR(w, w) + wR(w, w)) ( 1 — — 
2 \ \w\* 

k — w z R(w, w) + w z R(w, w). 



Expanding the integral quantities and the jump condi- 
tions into the asymptotic approximations: 



P = 
K = 



p(w, w)dz = P a + eP% + 0(e 2 ), 

) 

k(w, w)dz = eKi + e 2 K 2 + 0(e 3 ) 



and 



U?(s)[p]t = -4e 2 K (aia 2 +6ifo 2 ) + 0(e 3 ), 
Ue(s)[r}t = -2e 2 (a 1 6 1 +4 K 2 a 2 6 2 ) + 0(e 3 ) 



Ut(s)[r. 



= 0(e% 



we find the explicit expressions 
Pa = 2vn — 2 arctan ^— ^ , 

Uf(s)P 1 =2Kq-qd u P 
K x =A K {\- V 2 )U' e (s), 



Zvq-0 
2k 



d K P = 2U' e (s), 



U 2 e {s)K 2 = -8vKqU' e (s) - qd v K x + 
= -&v{U' e {s)) 2 



3isq — 
2k~ 



d R Ki 



and 



Ut(s)[-vp + r}± = -2e 2 v (U' e {s)) + 0(e 3 ). 



We note that the expressions for Pq and K\ are found by 
using wo (z) in Eq. I|25|l as a function of two independent 
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parameters k and v. This technical trick allows us to 
compute the correction terms Pi and K 2 by using partial 
derivatives of Pq and K\ in n and v. The final expressions 
use the relation k = y/1 — v 2 as well as the previously 
found relations |@5J|, (gHJ, and (JB3). 

Substituting all expansions into the balance equation 
((64(1 and integrating on z £ R, we obtain an extended 
main equation in the form: 



i> - n 2 U' c {s) + 



evU'J{s) 
2nU e {s) 



0(6 2 )- 



Using the relation s = vU e (s), the main equation can be 
rewritten for s(T): 



Ue(s)U' e ( S ) 



esUt'(s) 



0(e 2 ). 



(65) 



2/ct/ e (s) 

Furthermore, using the Tomas-Fermi approximation 
U e (s) = vl — s 2 , the extended dynamical equation is 
rewritten in the final form: 



2k(1 



+ 0(e 2 ), 



(66) 



where 



and s 2 + s 2 < 1. We note that the perturbation term Ri 
in the right-hand-side of the second-order problem 1(28(1 
does not contribute to the main equation (|66|l . since all 
associated integrals in the correction term Ki are zero 
due to symmetry in z € R. 

Second-order corrections 0(e 2 ) to the extended equa- 
tion (|66|l can be incorporated to the asymptotic theory 
if the second-order problem ((28(1 is solved explicitly and 
the third-order problem associated to the perturbed NLS 
equation 1(241) is analyzed. It is beyond the scope of our 
manuscript to derive the error bounds between the solu- 
tion of the perturbed NLS equation (|2~4fl and the trun- 
cated stationary solution of the first-order and second- 
order problems l(27|l and 1(28(1, 



VI. FAILURE OF THE FORMAL 
PERTURBATION THEORY 

We shall consider the perturbed GP equation in the 
form 1(2L)|) . where the ground state U e (x) is truncated by 
the Thomas-Fermi approximation (JTSJ, such that U e — 

VT^l 2 

explicit form: 



The perturbed GP equation 1(20(1 takes the 



iw t + -w xx + (1 - \w\ 2 )w = R(w, to), (67) 



where 



R(w, to) = e 2 a; 2 (l - \w\' z )w + 



1 — £ 2 2! 2 



to x . 



The GP equation 
R(w,w) = 0: 



l|67|l has the exact solution for 



w(x,t) = woirj) — fctanh(/c?7) +iv, 



where 



T) = X 



s(T) 



v(T) 



T = et. 



(68) 



(69) 



If Wo(r]) is a steadily traveling dark soliton, parameters 
are constant, where k = Vl — v 2 is the amplitude and 
v is the speed. We assume that the coordinate s(T) of 
the dark soliton changes adiabatically under the small 
perturbation R(w, w) ^ and show that a formal per- 
turbation theory fails to capture the correct dependence 
of the frequency of oscillations. Using the same renor- 
malized momentum equation (|64J) . we obtain from l|67(l 
the leading-order balance equation for the renormalized 
momentum of the dark soliton 0, Q, ^| : 



dP s 
' dT 



w' (x) (R + R) (too, wo)dx, (70) 



where 



P, = r 



(w w' - ww' ) ( 1 - 



I wo I 



dx 



2vk — 2arctan ( — 

v 



(71) 



such that 



dP s 
: dT 



4ek 2 k 
VT^k 2 ' 



The same equation occurs in the Lagrangian averaging 
technique applied to the perturbed GP equation 1)6 7|) (see 
[25| and references therein). The integrals in the right- 
hand-side of 1)70(1 can be evaluated at the leading-order 
approximation as e — ► 0, when the scaling 1)69(1 is used. 
As a result, we have 



< - / x 2 (l - \w \ 2 )(\w \ 2 )'dx = -ek 3 s + 0(e 3 ), 



2e 2 



rdx = — 



8 ek 3 t 



3 1 



+ 0(e 3 ). 



Using the leading-order approximation s = \/\ — k 2 , we 
close the main equation for s(T) as follows: 



(72) 



The main equation ((72(1 represents the adiabatic approx- 
imation for dynamics of dark solitons, in neglecting of 
the radiative effects (see |25|)- In the limit of black 
soliton, when s is small, the equation for an anharmonic 
oscillator ((72(1 approaches the equation for a harmonic os- 
cillator ((40(1 . However, for a dark soliton of arbitrary am- 
plitude, the anharmonic oscillator equation ((72|) is differ- 
ent from the harmonic oscillator equation 140|) , although 
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it represents the same asymptotic limit of adiabatic os- 
cillations of the dark soliton. 

The paradox above has a simple resolution. The first 
term in the right-hand-side of the perturbed GP equation 
(|67|l is not a small perturbation for the dark soliton i|68|) 
under the scaling H69(l in the limit e — *■ 0. Indeed, e 2 x 2 — 
s 2 — 2esi] + e 2 i] 2 , such that the correct perturbed GP 
equation in the variables l|69|) takes the form: 

iwt — ivwr) + —w vri + (1 — s 2 )(l — \w\ 2 )w = R(w, w), (73) 
where 

R = -2es V {l - \w\ 2 )w + -^w v + 0(e 2 ). 

1 — s £ 

However, the dark soliton w — Wo(f]) is n0 longer a so- 
lution of the left-hand-side of the GP equation (|73[1 and 
the renormalization of the variable r\ is required before 
formal application of the perturbation theory. The renor- 
malization of the GP equation is developed in the main 
part of this paper. We conclude that the formal appli- 
cation of the perturbation theory to the perturbed GP 
equation (|(j7|l fails to recover an accurate dependence of 
frequency of dark soliton oscillations from the amplitude 
of the dark soliton. 



VII. CONCLUSIONS 

In this paper, we have analyzed the oscillations of dark 
solitons in trapped atomic Bose-Einstein condensates. 
We have considered a repulsive quasi-one-dimensional 
condensate, described by a Gross-Pitaevskii equation 
with a parabolic external trapping potential. An asymp- 
totic multi-scale expansion method has been developed 
in the limit of the weakly trapped condensate. To 
the leading-order of approximation (i.e., neglecting the 
inhomogeneity- induced sound emission by the soliton), 



the soliton motion is harmonic, with a constant frequency 
depending on the trap strength. This result bridges ear- 
lier predictions obtained by different approaches in two 
limiting cases, namely the adiabatic perturbation theory 
for nearly black solitons [9j and the Korteweg-deVries 
approximation for shallow solitons 01 ■, an d is also in 
agreement with the recent prediction based on the semi- 
classical Landau dynamics of the dark soliton |T7| . On 
the other hand, to first-order approximation, the radi- 
ation (sound waves) emitted by the soliton due to the 
inhomogeneous background is taken into account. It is 
shown that radiation plays an important role in the soli- 
ton motion, as it responsible to an anti-damping effect, 
resulting in the increase of the amplitude of oscillations 
and decrease of the soliton amplitude. Energy loss of the 
soliton due to the emission of radiation is approximately 
found to follow an acceleration-square law, in agreement 
with previous numerical observations [l9( . We have also 
compared the results of the presented multi-scale expan- 
sion technique with the formal perturbation theory for 
dark solitons. We have found that a formal application 
of the perturbation theory fails to incorporate the cor- 
rect dependence of the frequency of the dark soliton os- 
cillation on the soliton amplitude (velocity). Numerical 
results have been found to be in good agreement to the 
analytical predictions, within the applicability intervals 
of the analytical assumptions. Finally, it should be no- 
ticed that the presented technique can also be used for 
the study of dark soliton dynamics in other relevant (in- 
homogeneous) systems, such as optical lattices and su- 
perlattices. 
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